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SUMMARY 


A new Navier-Stokes solver is developed by combining the 
efficiency of the LU-SSOR scheme and the accuracy of the flux- 
limited dissipation scheme. Application to laminar and turbu- 
lent flows and hypersonic flows proves the reliability of the 
new algorithm. 


1. INTRODUCTION 

An obvious way to accelerate convergence to a steady state 
is to increase the size of time step. Implicit schemes can 
take larger time steps than explicit schemes when the time 
step of an explicit scheme is restricted by stability rather 
than accuracy. Since the unfactored implicit scheme produces 
a large block banded matrix which is very costly to invert, 
the factored implicit schemes have been popular. However, 
using a large time step resulted in a large factorization 
error. Moreover, inherent stability problems occurred when 
the number of factors were greater than or equal to three. 
Nevertheless, several factored implicit schemes have proved to 
be useful [1 to 5] . 

An alternative way to get a steady-state solution is to 
solve the steady governing equations by the Newton method. 
Because of the rapid growth of the operation count with the 
number of mesh cells, the system was solved indirectly by the 
relaxation methods. Recently, a new relaxation method, the 
LU-SSOR scheme, was developed by combining the advantages of 
LU factorization and Gauss-Seidel relaxation. The vectorizable 
LU-SSOR scheme, which is based on central differences, requires 
scalar diagonal inversions [6]. 

It has long been recognized that upwind differencing can 
eliminate undesirable oscillations near shock waves. Stemming 
from the mathematical theory of scalar conservation laws, 

Harten proposed the concept of total variation diminishing 
(TVD) schemes [7]. TVD schemes preserve the monotonicity of 
an initially monotone profile, because the total variation 
would increase if the profile ceased to be monotone. Second 
order schemes can be constructed by using multipoint extrapo- 
lation formulas to estimate the numerical flux, or by adding 
higher order dissipative terms. In either case flux limiters 
are then needed to control the signs of the coefficients of a 
semi-discrete approximation to the hyperbolic system. However, 
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TVD schemes in two or more space dimensions are at most first 
order accurate. Nevertheless, two-dimensional results using 
one-dimensional second order TVD schemes and dimensional split- 
ting showed sharp resolution of discontinuities without oscil- 
lations. Jameson constructed an efficient scheme which can be 
converted to a TVD scheme in the case of a scalar conservation 
law using flux limited dissipation [8], 

In this paper, a new Navier-Stokes solver is developed by 
combing the efficiency of the LU-SSOR scheme and the accuracy 
of the flux limited dissipation scheme. Numerical examples 
include laminar and turbulent airfoils and a hypersonic inlet. 

2. THE NAVIER-STOKES EQUATIONS 

The Navier-Stokes equations represent gas flow in thermo- 
dynamic equilibrium. Let t, p, E, T, and p be time, den- 
sity, total energy, temperature, and pressure; u and v 
Cartesian velocity components; and x and y Cartesian coor- 
dinates. Then for a two-dimensional flow these equations can 
be written as 


3W 3F 3G 
3t + dx + dy 



3G 
v 

ay 


where W is the vector of dependent variables , and 
F and G are convective flux vectors: 


( 1 ) 


W = (p, pu, pv, pE)* \ 

F = [pu, pu 2 + p, pvu , u(pE + p)]*> (2) 

G = [pv, puv , pv 2 + p, v(pE + p)]*J 

Here * denotes the transpose of a matrix. The flux vectors 
for the viscous terms are 

/ 8T \ * 

F = [ 0, T , T ,UT +VT +kr~\ 

v l xx xy xx xy 3x I 

/ 3T \* 

G = 0, i , t ,ut +vt +k~) 

v y xy yy xy yy 3y / 

Here the viscous stresses are 


T 


XX 


= 2yu 

x 


T 

xy 


- f + V 

„ ( u y + v x ) 


and 


2 


2 

T = 2yv - ~ y(U + V ) 
yy y 3 ^ x y 

where y is the coefficient of viscosity and k is the 
coefficient of thermal conductivity. 


The pressure is obtained from the equation of state: 


P 


P(Y - 1) 



+ v 2 )| 


where y is the ratio of specific heats. 


( 3 ) 


3. THE LU-SSOR SCHEME 


A prototype implicit scheme for a system of nonlinear 
hyperbolic equations such as the Euler equations can be formu- 
lated as 

n+1 n „ /f n+l. ^ _ n+1. | 

W = W - (3At D^F(W ) + DyG(W ) 

- (l-(5) At jo x F(w n ) 

where D x and Dy are central difference operators that 
approximate d/3x and 3/3y. Here n denotes the time level. 
An enormous number of computations must be performed when the 
scheme is in this form because coupled nonlinear equations 
must be solved at each time step. Let the Jacobian matrices be 


D G(W 

y 




(4) 


A 

B 


3F ) 
dW f 

3G ( 

aw / 


(5) 


and let the correction be 

AW = W n+1 - W n 


The scheme can be linearized by setting 

F(W n+1 ) = F(W n ) + AAW + CKIIAWII 2 ) 

G(W n+1 ) = G(W n ) + BAW + 0( II AWII 2 ) 

where 0 is the order of the enclosed terms, and dropping 
terms of the second and higher orders to yield 


[1 + 3At(D x A + DyB) ] AW + At R = 0 (6) 
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where I is the identity matrix and R is the residual 


R = D x F(W n ) + D y G(W n ) 

If a constant 3 = 1/2, the scheme remains second-order accu- 
rate in time; for other values of P, the time accuracy drops 
to first order. 

The unfactored implicit scheme (Eq. (6)) produces a large 
block banded matrix that can be inverted only by performing a 
great many computations. In addition, a large amount of stor- 
age is required. If P = 1 the scheme reduces to a Newton 
iteration in the limit At ■* 

(D X A + DyB)6W + R = 0 (7) 

A diagonally dominant form of Eq. (7) is 

(D~A + + dV + D~B + + D + B~ ) 6W + R = 0 (8) 

x x y y 


where D and D are backward difference 

x y 

operators and D+ and D+ are forward difference 

operators. Here, two-point operators are used for 
steady flow calculations. A + f A“, B + f and B~ are constructed 
so that the eigenvalues of matrices are nonnegative and 

those of matrices are nonpositive: 


where 


A + = J(A + r A I) 
A’ - 1(A - r 4 I>, 
B* - J(B + r B I) 

B ' * J< B - r B IJ 

r > max ( |X I ) 
A - A 

r g > max ( | X g | ) 


(9) 


( 10 ) 


Here, X^ and X B represent eigenvalues of Jacobian 
matrices . 



Then, the LU-SSOR scheme for approximate Newton iteration 
becomes 


(D A + + D B + - A" - B )(D + A' + D + B 
x y x y 

+ A + + B + HW = - (r A + r B >(D x F + D y G) (11) 

which can be inverted in two steps as 

(D" A + + D - B + - A“- B") iW* = -(r. + r)(D F + D G) 
x y a b x y 

(D + A~ + D + B~ + A + + B + ) 6W = 6W* (12) 

x y 

For the Navier-Stokes equations, F and G are replaced 
by F - F v and G - G v in Eq. (11). That is, 


(dV + D y B + - A - B") (dV + dV ♦ k* + B + ) 

' - (r * + r B> { D X (F * V + °, <<S ' 

+ flux limited dissipation terms | 


&W 


(13) 


Since one-sided difference schemes are naturally dissipa- 
tive, no implicit smoothing is required on the left side. 

Only adaptive dissipation terms are explicitly added to the 
residual on the right side. It is interesting to note that 
the present numerical method eliminates the need for block 
diagonal inversions without using the diagonalization process. 
This is an especially desirable feature for the analysis of 
hypersonic reacting flows. The LU family of algorithms are 
fully vectorizable along i + j = constant lines on a vector 
computer. 

4. JAMESON’S TVD SCHEME 


Consider the scalar conservation law 


3u 

at 



= 0 


It is well known that the total variation 


TV . | |£| d* 

•'-00 


(14) 


(15) 
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i 


cannot increase. Suppose now that a multipoint semi-discrete 
approximation to equation (14) is expressed in the form 


du Q-l 


dt 


The discrete total variation is 


TV = E i u i - u i_i 


l=-00 

It can be shown that this will not increase if and only if 


and 


c_ x (i-l) > c _ 2 (i-2). . . > c Q (i-Q) > 0 
-c Q (i) > - c^(i+l) . . . > -c ^(i+Q-1) > 0 


Now consider the scheme 
du. 


dt + Ax (h i+l/2 h i-l/2 ) ° 


where h, , ^ is an approximation to the flux across the 
1+1/2 st th 

boundary between the (i+1) and i cells. Denoting f(u, 
by f^, define the numerical flux as 


h i+l/2 ~ 2 <f i+l + f i J + d i+l/2 


where d. . ^ is a dissipative flux. 
1 + 1/2 


d i+l/2 ” “i+1/2 




e i+3/2’ e i+l/2 > 2 e i+l/2 

+ B(e i.H/2- e * 


where e. . = p. , - p. for example. 
i+l/2 l+l i 


i-1/2^ 


B (p,q) is known as Roe’s min mod function which can be 
defined as 

B(p,q) = js(p) + s(q) | min (|P|,|q|) 

where 


s(p) = 


2 if P * 0 
if P < 0 


(16) 


(17) 


(18) 


) 


(19) 


(20) 


( 21 ) 
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Here 


( 22 ) 


“i+1/2 ■ min <1/2 ’ V k l“i*l/2 )<E ur V 

where k and k are the constants and 
o 1 


v l+l/2 max(v i+2’ v i+l ,v i ,v i-l ) 


V. 

1 


p. , 2P. P. n 

1+1 - 14 - 1-1 

P. - + 2P. + P. , 
l+l l l-l 


The spectral radius R can be estimated as the value of 
R = | Ayu - Axv | + cVax^ + Ay^ 


on the edge separating cells (i+l,j) and (i,j) in two- 
dimensions where c is the speed of sound. 


5. RESULTS 


The first test case was for viscous laminar flow past the 
NACA 0012 airfoil at Mach 0.5, Reynolds number 5000, and zero 
angle of attack. The adiabatic wall boundary condition was 
used at the body surface. Calculations were performed on a 
stretched 192 by 48 C-mesh. Figure 1 shows the Mach number 
contours while Fig. 2 shows velocity vectors. 

The next case was for viscous turbulent flow past the RAE 
2822 airfoil at Mach 0.73, Reynolds number 6.5 million, and 
2.79° angle of attack. The Reynolds - averaged Navier-Stokes 
equations were solved using a Baldwin-Lomax turbulence model. 
Transition was fixed at 3 percent chord. Mach number contours 
are shown in Fig. 3 where the dashed line denotes the sonic 
line. It is interesting to note that present TVD scheme pro- 
duces identical solutions to the adaptive dissipation results 
in Ref. 6 in the case of the weak shock waves. 

The last cases were for inviscid hypersonic flows past an 
inlet. Calculations were done on a uniform 54 by 32 H-mesh. 
Figures 4 and 5 show the Mach number contours for Mach 5 and 
10 flows respectively. When compared to the results of Ref. 9 
using the adaptive dissipation, present TVD scheme results 
show sharper resolution and significant improvement in both 
local and global accuracy. 
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Figure 1. - Mach number contours for viscous laminar flow. 
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FIGURE 2. - VELOCITY VECTORS FOR VISCOUS LAMINAR FLOW. 
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